# Lehmann and Masterson
# Aid and Anti-Refugee Violence
# Bandwidth Variation Plot


####
# Table of Contents -------------------------------------------------------
# 0. Preamble

# 1. Data cleaning

# 2. Analysis
#  2.1 Physical assault
#  2.1.1.	linear
#  2.1.2.	linear drop outlier bin

#  2.2 Verbal assault
#  2.2.1. linear
#  2.2.2. linear drop outlier bin

# 3. Plot



########################################################
# 0. Preamble -------------------

## Clear all objects
remove(list = ls())
gc()

# Uncomment to install packages
# install.packages(c("foreign", "sandwich",  "lmtest",  'tidyverse'))

packages <- c("foreign",  "sandwich", "lmtest", 'tidyverse')
lapply(packages, require, character.only = T)

data <- read.dta("C:/Users/chris/Dropbox/TRASH BIN/4may/LehmannMastersonAPSR2020_replication/LehmannMastersonAPSR2020_data.dta")

# Define cluster-robust standard error function
clx.fm <-
  function(fm, dfcw, cluster){
    library(sandwich)
    library(lmtest)
    M <- length(unique(cluster))
    N <- length(cluster)
    dfc <- (M/(M-1))*((N-1)/(N-fm$rank+M))  
    u  <- apply(estfun(fm),2,
                function(x) tapply(x, cluster, sum))
    vcovCL <- dfc*sandwich(fm, meat=crossprod(u)/N)*(fm$df / (fm$df - (M -1)))
    coeftest(fm, vcovCL) }


########################################################
# 1. Data cleaning

# Drop item non-response
data[(data[, 'q207'] == -96 | data[, 'q207'] == -97), 'q207'] <- NA
data[(data[, 'q208'] == -96 | data[, 'q208'] == -97), 'q208'] <- NA
data <- data %>% filter(is.na(data$q207)==0 & is.na(data$q208)==0) 

# Subset to respondent rows (c101==1)
data <- data %>% filter(data$c101 == 1) 

# Create the treatment dummy:
data$treat <- as.numeric(data$elevation_oct13>=500)

# Create polynomial functions:
data$dist <- data$elevation_oct13-500  # linear
data$dist_treat <- I(data$dist*data$treat)      # linear spline

# Create violence variables
# Create a 1/0 dummy 
data$violence <- (data$q208<3)  
data$verbal <- (data$q207<3)

# Rename baseline geographic area code
data <- data %>% rename(cas_code = cas_code_apr14)

# Define data without outlier bin
data_drop <- data %>% filter(dist < -10 | dist >= 0)

# Define differences for x axis
alt_seq <- seq(from = 0, to = 30, by = 1)

#####################################################
# 2. Analysis

######################################################
########### 2.1 Physical Violence Results
results_violence_linear <- data.frame(matrix(NA, nrow = length(alt_seq), ncol = 6))
results_violence_linear_drop <- data.frame(matrix(NA, nrow = length(alt_seq), ncol = 6))
colnames(results_violence_linear) <- c('alt_seq', 'beta', 'lwr90', 'upr90', 'lwr95', 'upr95')
colnames(results_violence_linear_drop) <- c('alt_seq', 'beta', 'lwr90', 'upr90', 'lwr95', 'upr95')

results_violence_linear[, 1] <- 50-alt_seq
results_violence_linear_drop[, 1] <- 50-alt_seq

#  2.1.1.	linear
for(i in 1:length(alt_seq)){   # each testing at i altitudes
  data_alt <- data %>% filter(abs(dist)<=50-alt_seq[i])
  
  fit_violence_pre_holder <- lm(violence ~ dist + treat + dist_treat, data = data_alt)
  fit_violence_holder <- clx.fm(fm = fit_violence_pre_holder,  cluster = data_alt$cas_code)

  lwr90_holder <- fit_violence_holder["treat", "Estimate"] - 1.645*fit_violence_holder["treat", "Std. Error"]
  upr90_holder <- fit_violence_holder["treat", "Estimate"] + 1.645*fit_violence_holder["treat", "Std. Error"]
  lwr95_holder <- fit_violence_holder["treat", "Estimate"] - 1.96*fit_violence_holder["treat", "Std. Error"]
  upr95_holder <- fit_violence_holder["treat", "Estimate"] + 1.96*fit_violence_holder["treat", "Std. Error"]
  
  results_violence_linear[i, 2:6] <-c(fit_violence_holder["treat", "Estimate"], lwr90_holder, upr90_holder, lwr95_holder, upr95_holder)
  
}

#  2.1.2.	linear -- drop outlier bin
for(i in 1:length(alt_seq)){   # each testing at i altitudes
  data_drop_alt <- data_drop %>% filter(abs(dist)<=50-alt_seq[i])
  
  fit_violence_pre_holder_drop <- lm(violence ~ dist + treat + dist_treat, data = data_drop_alt)
  fit_violence_holder_drop <- clx.fm(fm = fit_violence_pre_holder_drop,  cluster = data_drop_alt$cas_code)
  lwr90_holder_drop <- fit_violence_holder_drop["treat", "Estimate"] - 1.645*fit_violence_holder_drop["treat", "Std. Error"]
  upr90_holder_drop <- fit_violence_holder_drop["treat", "Estimate"] + 1.645*fit_violence_holder_drop["treat", "Std. Error"]
  lwr95_holder_drop <- fit_violence_holder_drop["treat", "Estimate"] - 1.96*fit_violence_holder_drop["treat", "Std. Error"]
  upr95_holder_drop <- fit_violence_holder_drop["treat", "Estimate"] + 1.96*fit_violence_holder_drop["treat", "Std. Error"]
  
  results_violence_linear_drop[i, 2:6] <-c(fit_violence_holder_drop["treat", "Estimate"], lwr90_holder_drop, upr90_holder_drop, lwr95_holder_drop, upr95_holder_drop)
  
}

######################################################
################## 2.2 Verbal Assault Results

results_verbal_linear <- data.frame(matrix(NA, nrow = length(alt_seq), ncol = 6))
results_verbal_linear_drop <- data.frame(matrix(NA, nrow = length(alt_seq), ncol = 6))
colnames(results_verbal_linear) <- c('alt_seq', 'beta', 'lwr90', 'upr90', 'lwr95', 'upr95')
colnames(results_verbal_linear_drop) <- c('alt_seq', 'beta', 'lwr90', 'upr90', 'lwr95', 'upr95')

results_verbal_linear[, 1] <- 50-alt_seq
results_verbal_linear_drop[, 1] <- 50-alt_seq

#  2.2.1. linear
for(i in 1:length(alt_seq)){   # each testing at i altitudes
  data_alt <- data %>% filter(abs(dist)<=50-alt_seq[i])
  
  fit_verbal_pre_holder <- lm(verbal ~ dist + treat + dist_treat, data = data_alt)
  fit_verbal_holder <- clx.fm(fm = fit_verbal_pre_holder,  cluster = data_alt$cas_code)
  lwr90_holder <- fit_verbal_holder["treat", "Estimate"] - 1.645*fit_verbal_holder["treat", "Std. Error"]
  upr90_holder <- fit_verbal_holder["treat", "Estimate"] + 1.645*fit_verbal_holder["treat", "Std. Error"]
  lwr95_holder <- fit_verbal_holder["treat", "Estimate"] - 1.96*fit_verbal_holder["treat", "Std. Error"]
  upr95_holder <- fit_verbal_holder["treat", "Estimate"] + 1.96*fit_verbal_holder["treat", "Std. Error"]
  
  
  results_verbal_linear[i, 2:6] <-c(fit_verbal_holder["treat", "Estimate"], lwr90_holder, upr90_holder, lwr95_holder, upr95_holder)
  
}

#  2.2.2. linear -- drop outlier bin
for(i in 1:length(alt_seq)){   # each testing at i altitudes
  data_drop_alt <- data_drop %>% filter(abs(dist)<=50-alt_seq[i])
  
  fit_verbal_pre_holder_drop <- lm(verbal ~ dist + treat + dist_treat, data = data_drop_alt)
  fit_verbal_holder_drop <- clx.fm(fm = fit_verbal_pre_holder_drop,  cluster = data_drop_alt$cas_code)
  lwr90_holder_drop <- fit_verbal_holder_drop["treat", "Estimate"] - 1.645*fit_verbal_holder_drop["treat", "Std. Error"]
  upr90_holder_drop <- fit_verbal_holder_drop["treat", "Estimate"] + 1.645*fit_verbal_holder_drop["treat", "Std. Error"]
  lwr95_holder_drop <- fit_verbal_holder_drop["treat", "Estimate"] - 1.96*fit_verbal_holder_drop["treat", "Std. Error"]
  upr95_holder_drop <- fit_verbal_holder_drop["treat", "Estimate"] + 1.96*fit_verbal_holder_drop["treat", "Std. Error"]
  
  results_verbal_linear_drop[i, 2:6] <-c(fit_verbal_holder_drop["treat", "Estimate"], lwr90_holder_drop, upr90_holder_drop, lwr95_holder_drop, upr95_holder_drop)
  
}


# 3. Plot
# NB: Change pdf destination as needed
pdf(file = "C:/Users/chris/Dropbox/TRASH BIN/4may/LehmannMastersonAPSR2020_replication/Figure2.pdf", width = 12, height = 9)

par(mfrow=c(2, 2)) 

######### Verbal  Full Data

plot(x = alt_seq, y = results_verbal_linear[, 'beta'], 
     ylim = c(-.475, .1),  
     type = 'l',
     ylab =  'Estimated Effect of Aid',
     xlab = 'Bandwidth (in meters)',
     cex.lab = 1.3,
     xaxt = 'n',
     main = "Probability of Verbal Assault (all obs)",
     cex.main = 1.4)
axis(1, at=c(0, 5, 10, 15, 20, 25, 30), labels=c(50, 45, 40, 35, 30, 25, 20))

polygon(x = c((alt_seq), rev((alt_seq))), 
        y = c(results_verbal_linear[, 'lwr95'], rev(results_verbal_linear[, 'upr95'])), 
        col = alpha('grey89', 0.5), border = NA)

polygon(x = c((alt_seq), rev((alt_seq))), 
        y = c(results_verbal_linear[, 'lwr90'], rev(results_verbal_linear[, 'upr90'])), 
        col = alpha('grey82', 0.5), border = NA)

lines(x = alt_seq, y = results_verbal_linear[, "beta"], 
      col = 'black', lwd = 1.5)

abline(h = 0, lty = 3, col = 'grey31')


######## Physical Linear Plot Full Data

plot(x = alt_seq, y = results_violence_linear[, 'beta'], 
     ylim = c(-.475, .1),  
     type = 'l',
     ylab = NA,
     xlab = 'Bandwidth (in meters)',
     cex.lab = 1.3,
     xaxt = 'n',
     main = "Probability of Physical Violence (all obs)",
     cex.main = 1.4)
axis(1, at=c(0, 5, 10, 15, 20, 25, 30), labels=c(50, 45, 40, 35, 30, 25, 20))

polygon(x = c((alt_seq), rev((alt_seq))), 
        y = c(results_violence_linear[, 'lwr95'], rev(results_violence_linear[, 'upr95'])), 
        col = alpha('grey89', 0.5), border = NA)

polygon(x = c((alt_seq), rev((alt_seq))), 
        y = c(results_violence_linear[, 'lwr90'], rev(results_violence_linear[, 'upr90'])), 
        col = alpha('grey82', 0.5), border = NA)

lines(x = alt_seq, y = results_violence_linear[, "beta"], 
      col = 'black', lwd = 1.5)

abline(h = 0, lty = 3, col = 'grey31')


##########################################################
######### Linear Plot -- Drop Outlier Bin
plot(x = alt_seq, y = results_verbal_linear_drop[, 'beta'], 
     ylim = c(-.475, .1), 
     type = 'l',
     ylab = 'Estimated Effect of Aid',
     xlab = 'Bandwidth (in meters)',
     cex.lab = 1.3,
     xaxt = 'n',
     main = "Probability of Verbal Assault (without outlier)",
     cex.main = 1.4)
axis(1, at=c(0, 5, 10, 15, 20, 25, 30), labels=c(50, 45, 40, 35, 30, 25, 20))

polygon(x = c((alt_seq), rev((alt_seq))), 
        y = c(results_verbal_linear_drop[, 'lwr95'], rev(results_verbal_linear_drop[, 'upr95'])), 
        col = alpha('grey89', 0.5), border = NA)

polygon(x = c((alt_seq), rev((alt_seq))), 
        y = c(results_verbal_linear_drop[, 'lwr90'], rev(results_verbal_linear_drop[, 'upr90'])), 
        col = alpha('grey82', 0.5), border = NA)

lines(x = alt_seq, y = results_verbal_linear_drop[, "beta"], 
      col = 'black', lwd = 1.5)

abline(h = 0, lty = 3, col = 'grey31')

###############################################################
######### Physical Linear Plot -- Drop Outlier Bin
plot(x = alt_seq, y = results_violence_linear_drop[, 'beta'], 
     ylim = c(-.475, .1), 
     type = 'l',
     ylab = NA,
     xlab = 'Bandwidth (in meters)',
     cex.lab = 1.3,
     xaxt = 'n',
     main = "Probability of Physical Violence (without outlier)",
     cex.main = 1.4)
axis(1, at=c(0, 5, 10, 15, 20, 25, 30), labels=c(50, 45, 40, 35, 30, 25, 20))

polygon(x = c((alt_seq), rev((alt_seq))), 
        y = c(results_violence_linear_drop[, 'lwr95'], rev(results_violence_linear_drop[, 'upr95'])), 
        col = alpha('grey89', 0.5), border = NA)

polygon(x = c((alt_seq), rev((alt_seq))), 
        y = c(results_violence_linear_drop[, 'lwr90'], rev(results_violence_linear_drop[, 'upr90'])), 
        col = alpha('grey82', 0.5), border = NA)

lines(x = alt_seq, y = results_violence_linear_drop[, "beta"], 
      col = 'black', lwd = 1.5)

abline(h = 0, lty = 3, col = 'grey31')

dev.off()

